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<^ ■ Abstract 

> 

^^ ■ The electronic structure of the superconducting surface sheath in a type- 

II superconductor in magnetic fields Hc2 < H < Hcs is calculated self- 
0^ . consistently using the Bogoliubov-de Gennes equations. We find that the 



pair potential A(x) exhibits pronounced Friedel oscillations near the surface, 



i-rt ■ in marked contrast with the results of Ginzburg-Landau theory. The local 

O ■ density of states near the surface shows a significant depletion near the Fermi 



energy due to the development of local superconducting order. We suggest 
that this structure could be unveiled by scanning-tunneling microscopy stud- 
ies performed near the edge of a superconducting sample. 
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The study of surface superconductivity was initiated some thirty years ago by the sem- 
inal work of Saint- James and de Gennes |1|], who predicted that a magnetic field parallel 
to a superconductor-vacuum interface would nucleate a superconducting "surface sheath" 
before the onset of superconductivity in the bulk of the material. By solving the linearized 
Ginzburg-Landau (GL) equations subject to the boundary condition that the normal deriva- 
tive of the order parameter vanish at the interface, they found an approximately Gaussian 
order parameter profile localized within a coherence length ^ of the interface; the critical 
field for this surface superconductivity is iJcs = 1-69 ifc2, where Hc2 is the bulk critical field 
for the Abrikosov flux-lattice phase |l|-0. This phenomenon has been confirmed by measure- 
ments which observe a vanishing surface resistance at fields above Hc2 (for a review of the 
early experiments, see Refs. [0]and |^). Due to the presence of the superconducting nucleus 
there is also a depletion of states near the Fermi energy, which is reflected in a suppression 
of the tunneling conductivity at low bias, an effect which has been observed in certain Pb-Bi 
and Sn-In alloys 0]. 

While the theory of surface superconductivity is quite complete within the framework 
of GL theory, there remain several interesting unanswered questions and problems which 
can only be addressed within a microscopic theory. (1) In a microscopic theory the nat- 
ural boundary condition for the quasiparticle wavefunctions is that they vanish at the 
superconductor-vacuum interface, so that superconducting pair potential also vanishes at 
the interface. If we naively apply this microscopic boundary condition to the macroscopic 
GL equations, and solve the linearized GL equations subject to the boundary condition that 
the order parameter vanish at the interface (rather than its derivative) we find H^s = Hc2', 
i.e., there is no surface superconductivity. However, we can not go back and infer that the 
microscopic equations do not possess surface superconducting solutions; by the same token, 
it cannot be taken for granted that the surface-superconducting solutions of GL theory (with 
the zero derivative boundary condition) prove the existence of surface-superconducting solu- 
tions of the microscopic theory. The two approaches, valid at different length scales, utilize 
different boundary conditions, and it is generally difficult to connect the two f^. Whether 



surface superconductivity exists within a realistic microscopic model of superconductors re- 
mains an open question. (2) The GL equations are derived using the quasiclassical phase 
approximation, which neglects the effects of Landau level quantization of the electronic 
states, and is valid in low fields. In very high magnetic fields the Landau level quantization 
can become important. Recent theoretical work has predicted de Haas- van Alphen oscil- 
lations in Hc2{T), as well as possible re-entrant superconductivity in very high magnetic 
fields P|. Might there be re-entrant surface superconductivity which precedes the re-entrant 
bulk superconductivity? Again, the answer would require a microscopic calculation which 
does not invoke the quasiclassical phase approximation. (3) On a parallel note. Landau level 
quantization in the presence of a surface produces magnetic edge states, which have been the 
subject of intensive study in the context of the integer and fractional quantum Hall effects 
0. Understanding the role that edge states play in surface superconductivity may help us 
in answering the basic question: Why is superconductivity favored at a surface? (4) In the 
mixed state of type-II superconductors the spatial variation of the pair potential can produce 
a rich structure in the local density of states (LDOS) @J^, which can be measured directly 
using a scanning-tunneling microscope (STM) []rU|. Analogous structure in the LDOS will be 
produced by the superconducting surface sheath, which should be observable using a STM. 
Such STM studies would provide the first direct image of the surface sheath, and provide 
valuable information about the pair potential profile and local electronic structure; all of 
the previous experimental studies have only explored averaged properties, such as the pair 
potential averaged over the sample [0-^. We note that previous attempts at a microscopic 
theory of surface superconductivity have been confined to analytical [jlT]-0 or numerical 



P^ solutions of the linearized gap equation. All of these works invoke the quasiclassical 
phase approximation, and make approximations which are equivalent to assuming that the 
derivative of the pair potential vanishes at the surface; they therefore do not address the 
issues which we have raised above. There have been no calculations of the LDOS. 

In this Letter we will discuss our first attempts at addressing some of the questions 
raised above by numerically solving the the microscopic Bogoliubov-de Gennes (BdG) equa- 



tions PI self-consistently, in a magnetic field with realistic boundary conditions at the 
superconductor-vacuum interface. We retain fully the Landau level (and edge state) struc- 
ture. To make the calculations computationally tractable we assume a two-dimensional ge- 
ometry, which is also a good approximation for many layered superconductors (see below). 
Our results show that the microscopic BdG equations do indeed admit a superconduct- 
ing solution localized near the surface, for fields H > Hc2- The pair potential vanishes at 
the surface, but rises rapidly and eventually looks like the GL solution; there is a narrow 
"boundary layer" near the surface in which the GL solutions break down. However, unlike 
the GL solution we find large amplitude Friedel-like oscillations in the pair potential. Our 
LDOS calculations show a suppression at low energies, in the regions where the pair poten- 
tial is a maximum; it should be possible to resolve this structure in an STM measurement. 
The remainder of this paper is organized as follows: Following a brief review of the BdG 
formalism, we will discuss our numerical methods. We will then present our results for the 
self-consistent pair potential and the LDOS, for a particular choice of parameters. Further 
details of these calculations will appear in Ref. |T^ . 



The Bogoliuhov-de Gennes Equations. The BdG equations for the quasiparticle am- 
plitudes Ui{r) and fi(r) with excitation energy e^ > (measured relative to the Fermi energy) 



are 
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with A(r) the pair potential. The single-particle electron Hamiltonian is 



K = ^ \-inV - -A(r)' 
c 



+ F(r) - Ef, (2) 
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where Ep is the Fermi energy, V{r) is the surface potential, and A(r) is the vector potential 
(we do not consider any effects of spin). Any effects of the band structure of the material are 
subsumed in the effective mass m*. The pair potential must be determined self-consistently 
from the solutions of the BdG equations, as 

A(r) = <? E <(r)«.(r)[l-2/(e,)], (3) 

ei<huiD 



where g is the BCS attractive couphng, ujd is the Debye frequency, and /(e) is the Fermi 
function. The vector potential must also be determined self-consistently using Ampere's 
law with the current density determined by the quasiparticle wavefunctions 0,^. Once 
the quasiparticle wavefunctions have been computed self-consistently, we can calculate the 
thermally broadened local density of states, 

N{r,E) = -J:[\u^{r)\'f{E-e,) 

i 

+ K(r)|V(i5 + e,)], (4) 

with /'(e) = Of /de. This quantity is proportional to the local differential tunneling conduc- 
tance which is measured in a STM experiment. 

We now take the magnetic field H = Hz to be parallel to the vacuum/superconductor 
interface at x = 0, with the superconductor occupying the half-space x > 0. As we will 
eventually be interested in modeling quasi-two dimensional materials, we will neglect dis- 
persion in the 2;-direction. Assuming the interface to be perfectly impenetrable, we have 
u{x = 0) = v{x = 0) = (and therefore A(a; = 0) = 0). In the Landau gauge A = (0, Hx, 0) 
the BdG equations are translationally invariant in the y-direction, and so we factor out this 
dependence as 

A(r) = e'2^°^/''A(a;), (5) 

Here /^ = hc/eH is the magnetic length, Xq the orbit center for the pair potential, — oo < 
Xo < c>o the orbit center for u (— xq is the orbit center for v), and n = 0, 1, 2, ... is a Landau- 
level index which counts the number of nodes of the wavefunctions; the sums in Eqs. (^ 
and (^) are over all (xo,n). We are left with a set of coupled one-dimensional equations 
for Uxo,n{x) and Vxo,n{x). Because of our boundary condition at x = the corresponding 
eigenvalues e„(xo,Xo) will depend upon the positions of the orbit centers, unlike the bulk 
case in which the energies are degenerate with respect to xq. The effort involved in solving 



these equations can be substantially reduced by finding both positive and negative energy 
solutions for Xq > and taking advantage of the transformation [^^ e -^ —e,u{r) -^ v*{r), 
f(r) —>■ —u*{r), to convert the negative energy solutions for a;o > into positive energy 
solutions for xq < 0. 

Method of solution. The BdG equations are solved iteratively, as follows. We start with 
an initial guess for the amplitude and the phase of the pair potential, taken from GL theory, 
for instance. We then fix the orbit center xq and calculate the wavefunctions and energies in 
the range < e < huo, by writing the BdG equations as a set of finite-difference equations. 
The lattice spacing is determined so that variations on the scale of Tc/kp can be resolved. 
The resulting matrix equations are sparse, and can be diagonahzed using standard packages 
(we use LAPACK). This process is then repeated for new values of Xq. The range of values of 
Xq is determined so that the highest energy states (of energy huij) are approaching their bulk 
behavior, i.e., becoming independent of xq, which occurs at xq = 1[2(Ef + hujD)/huJc\^^'^, 
with huc = heH/m*c the cyclotron energy. The spacing between these points is again 
determined by requiring that structure on the scale of Ti/kp can be resolved. Once all of 
the wavefunctions have been determined, the amplitude of the pair potential is recalculated 
from Eq. (^ by summing over xq and n. The phase of the pair potential is also recalculated 
by using the self-consistency condition for the vector potential p|,P]. The entire process is 
then repeated until the relative error in the order parameter between successive iterations 
is less than 0.02. 

Several cases were tested to determine the reliability of the algorithm. When A(r) =0 
we have reproduced the wavefunctions and spectrum for electrons in a constant magnetic 



field in the presence of an impenetrable surface |T6|. The eigenvalues for states with orbit 
centers at the surface {xq = 0, Xq = 0) were within 1% of those found analytically; for xq 
large the usual bulk Landau levels were obtained. We have also used several different initial 
guesses for the pair potential, including the GL form and a constant pair potential. The final 
results are insensitive the form of the initial pair potential, but convergence is expectedly 
slower for the constant pair potential. In the results shown here we have used a Gaussian 
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profile centered near the surface x = for the initial pair potential amplitude, and we have 
used for our initial Xq the value obtained in GL theory |0]. 

We have chosen to model a layered (i.e., two dimensional) material whose parameters 
obey the weak-coupling BCS relations. We took the Fermi surface to be cylindrical with 
"fT^x-v ^ ^^2 ^^^ ''^x-v ~ 2^^; where m is the electron mass. Assuming the zero temperature 
gap A(0) = 1.1 meV and the zero temperature coherence length ^(0) = 100 A yields 
a Fermi velocity of vp = 7rA(0)^(0)//i = 5.27 x 10^ cm/s, kp^ = 11.0 A, and Ep = 
15.8 meV. With huo = 15 meV and gN{0) = 0.30, the zero field critical temperature is Tc = 
7.25 K, and the zero temperature quasiclassical critical field is i/c2(0) = 0.722 0o/27r^^(O) = 
2.38 T [|1^. These parameters are similar to those used by Gygi and Schliiter in their study 
of the core structure of vortices in NbSe2 H], which showed good agreement with STM 



measurements ||T0| . Our Fermi energy is probably unrealistically low; it was chosen to make 
our computations tractable, and thus represents a compromise between numerical efficiency 
and realistic modeling. We do not expect any qualitative changes in our results for larger 
values of Ep. 

We have chosen to present results for a temperature of 2 K (T/Tc = 0.28), and a magnetic 
field of 4 T {H/H^2{T) = 1.83), so that / = 128 A and fiujc = 0.23 meV. With the parameters 
given above this means that we must keep a total of 134 bulk Landau levels. Using the criteria 
stated above, we need 350 lattice points in a sample 30/ wide. A second impenetrable wall 
is placed at x = 30/ to aid in normalization, but we ignore any nucleation at that surface. 
The finite difference version of the BdG equations is then an 700 x 700 matrix equation. 
The range of Xq is not limited to lie within the sample; indeed, the states with xq < are 
precisely the edge states. We use 800 orbit centers ranging from xq = —35/ to 35/, which 
ensures that all of the states in question have attained their bulk values. Convergence was 
reached after five iterations, requiring ten hours on an IBM RS 6000/370 workstation. 

Discussion of Results. Our result for the amplitude of the self-consistent pair potential 
is given in Fig. |^; this result is plotted against the GL theory result, taken from Ref. []TB[. 



The pair potential does indeed vanish at the material surface and exhibits large Friedel 



oscillations away from the surface, with a period which is approximately n/kp = 34 A. 
Such oscillations of the pair potential near a surface also occur in zero field JSUTEH, and in 
the vicinity of a magnetic impurity [^ , and are the result of pair-breaking by the surface 
and impurity. The self-consistent orbit center for the pair potential is Xq = 0.19/, which 
is close to the position of the first maximum of the amplitude of the pair potential. The 
existence of a surface sheath at these relatively high magnetic fields is not inconsistent with 
variational calculations at T = [|l^,|l^], which give Hcsi^O) > 1.95/7c2(0). We have repeated 
our calculations at a field of 4.5 T, and find that the maximum amplitude of the pair 
potential is decreased. Likewise, increasing the temperature results in a smaller amplitude. 
We have not been able to pin down Hc3{T) using our method, as achieving self-consistency 
becomes delicate when the pair potential is small. The phase boundary is best determined 
by direct numerical solution of the linearized gap equation |^ , which will also us to attack 
the question of re-entrant behavior in Hcs- 

The resulting change in the electronic structure can be seen in the thermally broadened 
local density of states, N{x,E), which we have plotted at constant energy. Fig. |^, and 
constant position. Fig. ^ In Fig. |] we see that at low energies the wavefunctions have 
a reduced amplitude in the vicinity of the maximum of the pair potential. Due to the 
presence of local superconducting order, in Fig. |^ we see that close to the surface there is a 
suppression in the density of states at low energies, with a corresponding enhancement at 
energies above the bulk gap. Farther from the surface we obtain the bulk normal density of 
states. We expect that this structure could be resolved by a STM measurement which would 
scan from the interior of a sample to a surface parallel to the applied magnetic field. Such 
a measurement would provide the first direct observation of the superconducting surface 
sheath. 
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FIGURES 
FIG. 1. The pair potential amplitude calculated from the BdG equations, compared with the 



GL result |18|, at H/Hc2{T) = 1.83 and T/T^ = 0.28. The BdG solution vanishes at the surface 



and exhibits strong Friedel oscillations as a result of the impenetrable wall at x = 0. 

FIG. 2. Thermally averaged local density of states for as a function of position (normalized to 
the normal local density of states), at H/ Hc2{T) = 1.83 and T /Tc = 0.28. For E = {} the electron 
states have been pushed away from the material surface by the nucleation of the pair potential. 
At higher energies {E = 1.1 meV) the effect becomes smaller, and vanishes altogether at very high 
energies. 

FIG. 3. Thermally averaged local density of states as a function of energy (normalized to 
the normal local density of states), at H/Hc2{T) = 1.83 and T/Tc = 0.28. Close to the surface 
(x = 0.1, 0.5) there is a significant depletion in the LDOS at low energies, as well as an enhancement 
at energies above the bulk gap of A(0) = 1.1 meV. When x = 1.5 the LDOS approaches its bulk 
normal state value. 
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